An explanation of the shape of the universal curve of the Scaling Law 
for the Earthquake Recurrence Time Distributions 
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This paper presents an explanation ol a possible mechanism underlying the shape of the universal 
curve of Scaling Law for Earthquake Recurrence Time Distributions. The presented simple stochastic 
cellular automaton model is reproducing the gamma distribution fit with the proper value of the 
parameter 7 characterizing Earth's seismicity and also imitates a deviation from the fit at the short 
interevent times, as observed in real data. 

Thus the model suggests an explanation of the universal pattern of rescaled Earthquake Recur- 
rence Time Distributions in terms of combinatorial rules for accumulation and abrupt release of 
seismic energy. 

PACS numbers: 91.30.Px, 91.30.Ab, 45.70.Ht, 02.50.Ga, 05.45.-a 

Keywords: recurrence times, universality, scaling law, stochastic cellular automata, avalanches, toy models 
of earthquakes, forest-fire models, Markov chains 



o- 
(Dp 

o: 
•I— I 



>: 

(N. 
^ ■ 

^ ■ 

(N. 



X: 



Analyzing seismic catalogs, Corral [H has determined 
that the probability densities of the waiting times be- 
tween earthquakes for different spatial areas and magni- 
tude ranges can be described by a unique universal dis- 
tribution if the time is rescaled with the mean rate of 
occurrence. 

To unify diverse observations the spatiotemporal anal- 
ysis was carried out as follows. Seismicity is considered as 
result of a dynamical process in which collective proper- 
ties are largely independent on the physics of the individ- 
ual earthquakes. Following Q, events are not separated 
into different kinds (foreshocks, mainshocks, aftershocks) 
nor the crust is divided into provinces with different tec- 
tonic properties. Then, a region of the Earth is selected, 
as well as temporal period and a minimum magnitude 
Mc (for conditions and other details see [J Events 
in this space-tinie-magnitude window are considered as 
a point process in time (disregarding the magnitude and 
the spatial degrees of freedom) and are characterized only 
by its occurrence time ti, with i — 1 . . . N{Mc)- Then the 
recurrence (or waiting) time is defined by Ti = ti—ti^i. 

The entire Earth has been analyzed by this method and 
it appears that different regions' probability densities of 
waiting times, rescaled by the mean seismic rate, as a 
function of the rescaled recurrence time collapse onto a 
single curve / [H: 



D(t-M,) = R{Mc)f{R{Mc)T), 



(1) 



where mean seismic rate R{Mc) is given by R{Mc) — 
N{Mc)/T, (here T is a total time into consideration), and 
recurrence-time probability density -D(r; Mc) is defined 
as D{t;Mc) = Pro6[r < recurrence time < r -I- drj/r. 
The so called scaling function / admits a fit in the form 
of a generalized gamma distribution 



ff'\0) = C6i'^-iexp(- 



(2) 



where 7 = 0.67 ± 0.005, /3 = 1.58 ± 0.15, S = 0.98 ± 0.05, 
C — 0.5 ± 0.1, and 9 — Rt is dimensionless recurrence 
time. The value of 6 can be approximated to 1. The 
present characterization of the stochastic spatiotempo- 
ral occurrence of earthquakes by means of a unique law 
would indicate the existence of universal mechanisms in 
the earthquake-generation process [1] . 

This paper presents an explanation of a possible mech- 
anism underlying the shape of the universal curve in 
terms of a cellular automaton model called Random 
Domino Automaton (RDA). The simple rules for evo- 
lution of the model, being a slowly driven system, rely 
on accumulation and abrupt release of energy only, which 
fit well to the described above procedure of neglecting in- 
dividual properties of earthquakes. We show that RDA 
reproduces the rescaled distribution of recurrence times. 

As can be seen from the original work I] as well as 
from further studies [3], results obtained from various 
earthquake catalogs show a deviation from the gamma 
distribution at the short interevent times. This holds 
from worldwide to local scales and for quite different tec- 
tonic environments. 

It is remarkable that the presented model reproduces 
also this deviation. Thus the model suggests an explana- 
tion of the universal pattern of rescaled Earthquake Re- 
currence Time Distributions in terms of its combinatorial 
rules for accumulation and release of seismic energy. 

So far, some insight of the origin of the gamma dis- 
tribution as well as examination the recurrence statistics 
of a range of cellular automaton earthquake models are 
presented in It is shown there, that only one model, 
the Olami-Feder-Christensen automaton, has recurrence 
statistics consistent with regional seismicity for a certain 
range of conservation parameter of that model. 

The Random Domino Automaton (RDA) was intro- 
duced as a toy model of earthquakes but can be 
also regarded as an extension of well known 1-D forest- 
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TABLE I: States of RDA for the size of the lattice iV = 5. 



state label 


example 




multiplicity 
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^ — ^ II 
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^ 1 1 1 • 
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3 


M> 1 1 1 1 • 1 • 


1 ^ 
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M> 1 1 1 • 1 1 • 


1 ^ 


5 


5 


^1 1 1 • 1 • 1 • 




5 


6 


^1 1*1 1 • 1 • 




5 


7 


{•{•{•• 




5 


8 






1 



fire model proposed by Drossel and Schwabl [10|. As a 
field of application of RDA we have already studied its 
relation to Ito equation Il|-|l3| and to integer sequences 
We point out also other cellular automata models 



17| giving an insight into diverse specific aspects of 
seismicity, including predictions. 

The RDA is characterized as follows: 

- space is 1-dimensional and consists of N cells; periodic 
boundary conditions are assumed; 

- cell may be empty or occupied by a single ball; 

- time is discrete and in each time step an incoming ball 
hits one arbitrarily chosen cell (the same probability for 
each one). The balls are interpreted as energy portions. 

The state of the automaton changes according to the 
following rule: 

• if the chosen cell is empty it becomes occupied with 
probability i/; with probability (1 — i^) the incoming ball 
is rebounded and the state remains unchanged; 

• if the chosen cell is occupied, the incoming ball provokes 
an avalanche with probability fi (it removes balls from hit 



cell and from all adjacent cells); with probability (1 — /i) 
the incoming ball is rebounded and the state remains 
unchanged. 

The parameter v is assumed to be constant but the 
parameter fx = is allowed to be a function of size i of 
the hit cluster. This extension with respect to Drossel- 
Schwabl model leads to substantial novel properties of 
the automaton. Note, only ratio of these parameters Hi/v 
affects properties of the automaton - changing of ji and 
v proportionally corresponds to a rescaling of time unit. 

The remarkable feature of the RDA is the explicit one- 
to-one relation between details of the dynamical rules 
of the automaton (represented by rebound parameters 
[lijv) and the produced stationary distribution of clus- 
ters of size i, which implies distribution of avalanches wi. 
It shows how to reconstruct details of the " macroscopic" 
behavior of the system from simple rules of " microscopic" 
dynamics. 

Various sizes iV of RDA can be considered in order to 
explain the shape of the universal curve of Scaling Law. 
It appears results for quite a small size N — 5 are enough 
to explain the idea and allow to keep the reasoning con- 
cise and detailed. RDA for a bigger size of the lattice 
behaves similar-like and the overall picture remains the 
same, as results from explanations given below. 

RDA is also a Markov chain 9]. It is convenient to 
define states i up to translational equivalence. Thus, in 
example, for = 5, instead of 2^, there are 8 states 
only - see Table HI Such space of states is irreducible, 
aperiodic and recurrent. Transition matrix p, where 
[p] ij = probability of transition i — > j, for = 5 is of 
the form 
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Stationary distribution is given by 
ti • p = w. 



(4) 



The evolution of the system is represented in Figure [T] 
Arrows between states i and j, with respective weights 
Wy , indicate possible transitions. A symbol L(j) depict 
an avalanche to state j. The density of the system is 
growing from left side (state 1 has density p = 0) to 
right side (up to density p = 1 for state 8). 



The expected time between two consecutive avalanches 
Tav may be expressed by various formulas [9] . For exam- 
ple 



(w) + 1 
1 - P ' 



(5) 



where {w) is the average avalanche size and Pr is the 
probability that the incoming ball is rebounded both 
form empty or occupied cell . 
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FIG. 1: A state diagram for RDA of size N = 5. Arrows with 
respective weights indicate possible transitions; those with 
avalanches are ended with symbol "L". A state is boxed, if it 
is possible to get it directly after an avalanche. 



The probabilities Vi of states i obtained from condition 
^ allow to determine the distribution of frequency fi of 
avalanche of size i, if rebound parameters fii/v are given. 
There exists also a systematic procedure of obtaining ap- 
proximate values of rebound parameters /Xi/iA which pro- 
duce requested distribution of avalanches [l8|. The ap- 
proximation comes from nonexistence of exact equations 
for (stationary) distribution of clusters Ui for sizes bigger 
than 4 (see [9j). We have used this procedure to obtain 
values ^J-i/v that give noncumulative inverse-power dis- 
tribution of avalanches presented in Table |lll 

To calculate the distribution of waiting times, each 
path starting from a state reached after an avalanche 
and ending with an avalanche is considered. There are 
42 such paths for size N = 5. Each path is assigned its 
respective probabilities, that a total passage time is equal 
to 1,2,... time steps. 

Respective weights Sk describing how often the system 
starts from initial state k are given by 



Sk 



J2i>k ^iPik 



For = 5 initial states are 1, 2 and 3. 
The expected time of stay in a state k is 



^av {k) 



- Pkk) - i = {I - PkkY 



(6) 



(7) 



Probability of stay in given state k for a time equal to 
i time steps is given by 



(i-l) 
^Pkk 



(1 - Pkk), 



(8) 



and all possible values are aggregated in a vector T*^ with 
i-th component equal to T/^. For path through two con- 
secutive states k and I, respective probability of time of 



TABLE II: Approximate values of rebound parameters fii and 
respective avalanche distribution Wi The parameter u — 0.25. 



i 


1 


2 


3 


4 


5 


Wi 


0.999060 
0.413247 


0.388232 
0.102851 


0.284504 
0.042587 


0.097650 
0.022351 


0.045810 
0.014306 



stay in both of them equal to j time steps is defined by 



rpkl ^ (rj^k ^ rj^ly^ ^ ^ yfey^I ^ 



(9) 



For a path through three states k,l,m we have T*^'™ = 
^rpk ^ yi) ^ and so on for longer paths. 

The probability rates Wij for transition i j where 
i ^ j are just normalized probabilities pij, namely 

= (10) 

l^j^i Pi-3 

,ik-i,ik there is assigned a 



Thus for a path Zi,j2, 
combined weight 



W 



llZ2...ifc-llfc 



(11) 



(12) 



as well as combined weigted time vector 

The ith component of the vector f7'i*2 •»fc gives a con- 
tribution to waiting time equal to i coming from a path 
ii, Z2, . • . , ik- Summing up those vectors for all possible 
paths we end with a distribution of waiting times. One 
can obtain also a distribution related to avalanches of 
chosen size. For example, if such sum is made for paths 
related to avalanches of size 2,3,4 and 5 only, a distribu- 
tion of waiting times related to avalanches of size bigger 
then 1 is obtained. 

Rebound parameters presented in Table HI] were cho- 
sen in order to obtain noncumulative distribution of 
avalanches in the form consistent with Gutenberg- Richter 
law. The exact value of power (here 2.1) does not affect 
results of the construction substantially. 

The system has average density p = 0.273885, average 
avalanche size (w) = 1.52458 and average time between 
avalanches Tav = 21.2027. The parameter Pr = 0.880932 
shows, that most of incoming balls are rebounded. Ex- 
pected times of staying in all states are presented in Ta- 
ble Hill Great majority of avalanches leads to empty state 
{Si — 0.755449), roughly every fifth avalanche leads to 
state 2 (52 = 0.205253), and roughly every twenty fifth 
to state 3 (^3 = 0.0392984). 

Figure [2] presents obtained distributions of waiting 
times up to 300 time steps. The upper curve is for 
avalanches of all sizes, the next is for avalanches of 
size bigger then 1, and so on. The lowest curve count 
avalanches of size 5 only. 
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TABLE III: Approximate expected stay times in states for 
the size of the lattice N = 5. 



state 


1 


2 


3 


4 


5 


6 


7 


8 




4.0 


2.5 


3.3 


1.8 


3.7 


2.2 


7.8 


21.8 



p 




1 10 100 

interevent time [time steps] 



smallest avalanches (also not recorded in real data). 

Note, that due to the incompleteness of the seismic 
catalogs in the short-time scale, usually real data are not 
displayed on plots for very short time intervals. Thus, 
the obtained theoretical curve, shown in Figure [U may 
be similarly cut for small times. If it is done for time, 
say, smaller then 10, it reflects the shape of real data. 

Thus, the presented model suggests that the origin of 
a universal curve is of combinatorial nature of accumula- 
tion and abrupt release of energy according to the rules 
depending on some parameters defining probabilities de- 
pendent on size of energy portions, as described above. 
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FIG. 2: A plot of distributions of interevent times for RDA 
with N=5. The upper line is for avalanches of all sizes, the 
curve below is for avalanches of size bigger then 1, the next 
is for sizes bigger then 2, and so on. The lowest curve counts 
only avalanches of size 5. The dashed line is a plot of fitted 
gamma distribution aj/'^~^'e~T; the solid line below is a plot 
of fitted exponential curve a'e~"S^ . 

Dashed line is a fitted gamma distribution ay'^^^-'e^'^ , 
where 7 = 0.67, b = 22.4 and c = 0.011. This fit 
is done for points with time coordinate from 60 (x^ = 
5.5346 ■ 10~^^). Values of the parameters b and c can be 
rescaled, depending on their relation to physical quan- 
tities (time, number of earthquakes). The parameter 7 
is a fixed parameter, with exactly the same value which 
characterize Earth's seismicity. The solid line below is a 
plot of fitted exponential curve a'e~^. 

Thus, the exponential part of the universal curve comes 
from distributions of biggest avalanches. In the presented 
example the biggest tav is for the state 8 containing single 
cluster of size 5 (see Table IIIip . Thus its contribution to 
the overall waiting time distribution dominates for big- 
ger times (compare formulas ([7]) and ([5])). Also state 7 
containing single cluster of size 4 contribute, but it is 
decaying more rapidly. 

The other part of the universal curve, comes from con- 
tributions of avalanches of smaller sizes. Its shape is a 
result of composition of many possible paths of the evo- 
lution, depicted in Figure [TJ For bigger sizes N there are 
much more possible paths (i.e. 1554 for N = 7) through 
states containing many clusters with comparable times 
tav Their composition flatten the curve. Moreover, cal- 
culation shows this effect produces a surplus (comparing 
to the gamma fit) for small waiting times, which is ev- 
ident in real earthquakes data [l|, Ij]- The size of the 
surplus can be reduced by omitting of a contribution of 
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